

if (strpos("`c(pwd)'", "xfan") > 0) {   // monarch
   local lvix = strpos("`c(pwd)'", "BP_MPI")
   local lvss = substr("`c(pwd)'", `lvix', .)
   local lvix = strpos("`lvss'", "/")
   global cwd_indiv = "/home/xfan/scratch/BP_output/"+substr("`lvss'", 1, `lvix'-1)+"/output"
}
else if (strpos("`c(pwd)'", "xfa561") > 0) {    // NCI
   local lvix = strpos("`c(pwd)'", "BP_MPI")
   local lvss = substr("`c(pwd)'", `lvix', .)
   local lvix = strpos("`lvss'", "/")
   global cwd_indiv = "/short/kw0/xfa561/BP_output/"+substr("`lvss'", 1, `lvix'-1)+"/output"
}
else if (strpos("`c(pwd)'", "/z3445921") > 0) {   // katana
   local lvix = strpos("`c(pwd)'", "uBP_MPI")
   local lvss = substr("`c(pwd)'", `lvix', .)
   local lvix = strpos("`lvss'", "/")
   global cwd_indiv = "/srv/scratch/z3445921/BP_output/"+substr("`lvss'", 1, `lvix'-1)+"/output"
}
else {
   global cwd_indiv = "."
}

if (strpos("`c(pwd)'", "BP_MPI_col")>0) {
   global gIcol = 1
}
else {
   global gIcol = 0 
}

global gvdata0 = 22 + ${gIcol}*4
*global gvdata1 = 50 
global gvdata_diff = 50
global gvdata10 = 65
global gvdata_ss = 60
global gvdata10_ss = 70 
global gvtNT = 76
global gvNT = 80

*global gvt_elas = ${gvdata10}
global gvt_elas = 70

global gvNT_incLE = 90



global alt_exo "exo" 
global alt_exo_id = 8
global alt_exo_des "Exogenous"
global alt_lbd "lbd"
global alt_lbd_id = 9
global alt_lbd_des "Learning-by-doing"

global alt_nodepatwork "nodepatwork"   // gcI_robust = 22
global alt_nodepatwork_id = 22
global alt_nodepatwork_des "No depreciation at work"

global gv_alt1 ${alt_exo}
global gv_alt1_id = ${alt_exo_id}
global gv_alt1_des ${alt_exo_des}
global gv_alt2 ${alt_lbd}
global gv_alt2_id = ${alt_lbd_id}
global gv_alt2_des ${alt_lbd_des}


*global alt_health "health"
*global alt_health_id = 20
*global alt_pt "pt"
*global alt_pt_id = 30

global alt_pthealth "pthealth"
global alt_pthealth_id = 40


global gvparnum = 25
global gvparnum_pthealth = ${gvparnum} + 8 


if (strpos("`c(pwd)'", "exo") > 0) {
   global alt_name "_exo"
}
else if (strpos("`c(pwd)'", "lbd") > 0) {
   global alt_name "_lbd"
}
else if (strpos("`c(pwd)'", "nodepatwork") > 0) {
   global alt_name "_nodepatwork"
}
else if (strpos("`c(pwd)'", "pthealth") > 0) {
   global alt_name "_pthealth"
}
else {
   global alt_name ""
}

* define the program to infile moments
cap program drop prog_MSM_moments_infile
program define prog_MSM_moments_infile
   local lvfortranmoments lfpr_f lnw_f sdlnw_f lnw_fd_f alfpr1to0_f alfpr0to1_f C_f ssa_f

   infile `2' t `lvfortranmoments' lfpr_diff_E2G_f lfpr_diff_G2B_f lfpr_diff_B2D_f lfpr_pt_f plnw41_f `3' using "`1'", clear
   sum t, detail
   local lvtmin = r(min)
   local lvtmax = r(max)
   if (`lvtmin'<10 | `lvtmax'>100) {
       infile `2' t `lvfortranmoments' lfpr_diff_E2G_f lfpr_diff_G2B_f lfpr_pt_f plnw41_f `3' using "`1'", clear
       gen lfpr_diff_B2D_f = .
   }

   foreach iv in `lvfortranmoments' {
      replace `iv' = . if `iv'<0
   }

   replace lnw_fd_f = . if lnw_fd_f<=0
   replace lnw_f = . if lnw_f<=0

   *replace ssa_f = . if t<${gvdata_ss}
   foreach iv in lfpr_diff_E2G_f lfpr_diff_G2B_f lfpr_diff_B2D_f {
      replace `iv' = . if t<${gvdata_diff}
      replace `iv' = . if `iv'<-990
   }

   foreach iv in lfpr_pt_f {
      replace `iv' = . if `iv'<-990
   }

   compress

   sort t
   label var t "Age"

   *keep if t>=${gvdata0} & t<=${gvdata10}
end


cap program drop prog_MSM_indiv_infile
program define prog_MSM_indiv_infile
   quiet infile id t A H AIME C I L labor EV LogitError XIerror HEALTHerror V0 V1 probV0 FS_Mt_err FS_spinc s_irecss irecss ihealth iFS ihet obs_start obs_end b_output b_sstax b_medicaretax b_finctax b_ssb b_ssdiinc Vp probVp using `1', clear
   compress
   count if V0>=.
   count if V1>=.
   gen olabor = 1 - L
   replace olabor = . if olabor <=0
   gen inc = H*labor
   gen lnw = log(H*labor/olabor) if olabor>0 & olabor<. & labor>0
   replace lnw = . if lnw<=-5.7
    
   gen lfpr_pt = (L<0.9 & L>0.2)
   gen lfpr_ft = (L<=0.2)

   *gen work = lfpr if lfpr_pt==0  // work does not include pt, same as in SIPP
   gen lfpr = (L<0.9) if lfpr_pt == 0    // exclude part time 

   gen I_cond = I if L<0.9
   gen lnH = log(H)
   compress

   * social security application age
   sort id t
   by id: gen ssa = (irecss==1 & s_irecss==0)
   label var ssa "Social Security Application"

   * obtain wage with fe
   *sort id t
   *xi: xtreg lnw i.t if t>=${gvdata0} & t<=${gvdata10}, fe i(id)
   *predict lnw_fe_full if e(sample)
   *xi: xtreg lnw i.t if t>=obs_start & t<=obs_end, fe i(id)
   *predict lnw_fe if e(sample)

   xtset id t
   sort id t
   gen lnw_ft = lnw if lfpr_ft==1
   by id: gen dlnw_full = d.lnw_ft   // _full meaning the full life cycle
   gen dlnw = dlnw_full if t>=obs_start & t<=obs_end   
   // dlnw is focusing on the full time lnw only, as what we did in the data
   
   by id: gen dlnw_base = lnw_ft if !missing(dlnw[_n+1])
   by id: gen dlnw_full_base = lnw_ft if !missing(dlnw_full[_n+1])
   drop lnw_ft

   
   *xi: xtreg lnw i.t if t>=${gvdata0} & t<=70, fe i(id)
   *predict lnw_fe if t>=${gvdata0} & t<=70
   *drop _It_??

   * E to U, U to E
   gen lfpr_next = .
   gen lnw_next = .
   *gen work_next = .
   sort id t
   by id: replace lfpr_next = lfpr[_n+1] if t[_n+1] == t+1
   by id: replace lnw_next = lnw[_n+1] if t[_n+1] == t+1
   *by id: replace work_next = work[_n+1] if t[_n+1] == t+1


   * adult equivelent consumption
   gen C_orig = C
   replace C = C / 1.34 if iFS==2 | iFS==3

   * the human capital of workers;
   gen H_worker = H if olabor>0 & olabor<. & labor>0
   gen lnH_worker = log(H_worker) if !missing(H_worker)
   *replace lnH_worker = . if lnH_worker<=-5.7
  
   compress

   sort id t
  
end


* for plot

global gytitle_lfpr "Employment Rates (ER)"
global gpos_lfpr = 6

global gytitle_lnw "Mean Log Measured Wages"
global gpos_lnw = 6

*global gytitle_lnw_fe "Mean Log Measured Wages (FE)"
*global gpos_lnw_fe = 6
global gytitle_lnw_fd "Mean Log Measured Wages (FD)"
global gpos_lnw_fd = 6


global gytitle_sdlnw "S.D. Log Measured Wages"
global gpos_sdlnw = 6

global gytitle_C "Adult Equivalent Consumption"
global gpos_C = 6

if (${gIcol}==0) {
    global gyrange_lfpr "0.2 1"
    global gylabel_lfpr "0.2(0.2)1"

    global gyrange_lnw "2.3 2.8"
    global gylabel_lnw "2.3(0.1)2.8"

    *global gyrange_lnw_fe "2.3 2.8"
    *global gylabel_lnw_fe "2.3(0.1)2.8"
    global gyrange_lnw_fd "2.3 2.8"
    global gylabel_lnw_fd "2.3(0.1)2.8"

    global gyrange_C "6 12"
    global gylabel_C "6(2)12"
}
else {
    global gyrange_lfpr "0.4 1"
    global gylabel_lfpr "0.4(0.2)1"

    global gyrange_lnw "2.6 3.4"
    global gylabel_lnw "2.6(0.2)3.4"

    *global gyrange_lnw_fe "2.6 3.4"
    *global gylabel_lnw_fe "2.6(0.2)3.4"
    global gyrange_lnw_fd "2.6 3.4"
    global gylabel_lnw_fd "2.6(0.2)3.4"

    global gyrange_C "8 18"
    global gylabel_C "8(2)18"
}

global gyrange_sdlnw "0 1"
global gylabel_sdlnw "0(0.25)1"

global gytitle_ssa "Social Security Application"
global gyrange_ssa "0 0.5"
global gylabel_ssa "0(0.1)0.5"
global gpos_ssa = 1


global gytitle_plnw41 "Wage Depreciation After Unemployment"
global gyrange_plnw41 "-0.12 0"
global gylabel_plnw41 "-0.12(0.03)0"
global gpos_plnw41 = 6

global glfprdiffupper 0.8
global gytitle_lfpr_diff "Difference in Employment Rates"
global gyrange_lfpr_diff "0 ${glfprdiffupper}"
global gylabel_lfpr_diff "0(0.2)${glfprdiffupper}"
global gpos_lfpr_diff = 6
global gytitle_lfpr_diff_E2G "Difference in ER: Excellent - Good"
global gyrange_lfpr_diff_E2G "0 ${glfprdiffupper}"
global gylabel_lfpr_diff_E2G "0(0.2)${glfprdiffupper}"
global gpos_lfpr_diff_E2G = 6
global gytitle_lfpr_diff_G2B "Difference in ER: Good - Bad"
global gyrange_lfpr_diff_G2B "0 ${glfprdiffupper}"
global gylabel_lfpr_diff_G2B "0(0.2)${glfprdiffupper}"
global gpos_lfpr_diff_G2B = 6
global gytitle_lfpr_diff_B2D "Difference in ER: Bad - Disabled"
global gyrange_lfpr_diff_B2D "0 ${glfprdiffupper}"
global gylabel_lfpr_diff_B2D "0(0.2)${glfprdiffupper}"
global gpos_lfpr_diff_B2D = 6

global gytitle_lfpr_pt "Employment Rates (ER) Part Time"
global gyrange_lfpr_pt "0 0.2"
global gylabel_lfpr_pt "0(0.05)0.2"
global gpos_lfpr_pt = 12 

global gytitle_lfpr_ft "Employment Rates (ER) Full Time"
global gyrange_lfpr_ft "0 1.0"
global gylabel_lfpr_ft "0(0.25)1.0"
global gpos_lfpr_ft = 12

cap program drop prog_parSSAx1000
program define prog_parSSAx1000
   local lvi = 1
   local lvx ``lvi''
   while ("`lvx'"!="") {
      replace `lvx' = `lvx' * 1000 if incode>=23 & incode<=25
      local lvin = `lvi' + 1
      local lvi = `lvin'
      local lvx ``lvi''
   }
end


* label CF
local lvCF0 baseline
local lvCF1 NoSSTest
local lvCF2 DelayNRAonly  //2--Delaying Social Security Normal Retirement Age (NRA)
local lvCF3 NoSSsystem
local lvCF4 NoUncertainty
local lvCF5 ReduceSSbenefit
local lvCF6 NoSSbenefit
local lvCF7 NoSStax
local lvCF8 IncreaseTax
local lvCF9 DelayNRAnUtil

local lvCF10 NoSS_wSSDI
local lvCF11 ReduceSSb_wSSDI
local lvCF12 NoSSb_wSSDI

local lvCF20 VaryPIA061r
local lvCF21 VaryPIA01r
local lvCF22 VaryPIA03r
local lvCF23 VaryPIA05r
local lvCF24 VaryPIA07r
local lvCF25 VaryPIA09r
local lvCF26 VaryAIMElast
local lvCF27 VaryAIMElastMax

        
